clean baseline

After the filtering steps and self-assignment results from 02- coupled with some investigation into mis-identified samples in 04-, it’s clear there are some fish that should be excluded (helvomaculatus, brevispinis, simulator) and others with incorrect meta data (maliger).

Here I clean up those things and see what the assignment looks like.

I also want to add some meta data to the samples so that I can understand any oddities - like with the hopkinsi…

library(tidyverse)
library(rubias)
# data
two_col <- read_csv("csv_outputs/sebastes_spp_baseline_two_col.csv")

Parsed with column specification:
cols(
  .default = col_integer(),
  sample_type = col_character(),
  repunit = col_character(),
  collection = col_character(),
  indiv = col_character()
)
See spec(...) for full column specifications.

# previous mis-assignments at the 95% level
mis <- read_csv("csv_outputs/baseline_misassignments95.csv")

Parsed with column specification:
cols(
  indiv = col_character(),
  collection = col_character(),
  repunit = col_character(),
  inferred_collection = col_character(),
  inferred_repunit = col_character(),
  scaled_likelihood = col_double(),
  log_likelihood = col_double(),
  z_score = col_double(),
  n_non_miss_loci = col_integer(),
  n_miss_loci = col_integer()
)

mis

# metadata
meta <- readRDS("../new_baseline_data/processed/meta-data-tibble.rds")
# from the nsf data, I know there are duplicates in the meta data

From this, swap the identity of the maliger (actually caurinus); remove the helvomaculatus, simulator, and brevispinis.

new2col <- two_col %>%
  mutate(repunit = ifelse(indiv == "R035252", "caurinus", repunit)) %>%
  mutate(collection = ifelse(indiv == "R035252", "caurinus", collection)) %>%
  filter(!repunit %in% c("helvomaculatus", "simulator", "brevispinis"))
new2col %>%
  write_csv("csv_outputs/new2col.csv")

Now try self-assignment to see where the rest of the misassignments fall.

# perform self-assignment of reference samples
new_self <- self_assign(new2col, gen_start_col = 5)

# correct assigments at the 95% threshold?
new_self %>%
  filter(repunit == inferred_repunit) %>%
  filter(scaled_likelihood > 0.95)

That is 1,579 of 1,647

1579/1647

[1] 0.9587128

95.9% accurate assignment at the 95% threshold. Great.

# and take a quick look at the misassignments
new_miss <- new_self %>%
  filter(inferred_repunit != repunit) %>%
  filter(scaled_likelihood > 0.95) # 9 high confidence mis-assignments
new_miss %>%
  filter(repunit %in% c("carnatus", "chrysomelas")) # 8 of those assignments are gopher/black-and-yellow

8/9

[1] 0.8888889

And for misassignments at a lower threshold?

new_self %>%
  filter(scaled_likelihood > 0.5 & scaled_likelihood < 0.95) %>%
  filter(!collection %in% c("carnatus", "chrysomelas"))

Warning in `[<-.data.frame`(`*tmp*`, is_list, value = list(`11` = "<>")) :
  replacement element 1 has 1 row to replace 0 rows

All 59 intermediate misassignments are gopher/black-and-yellow.

Misassignments without gopher/black-and-yellow

new_self %>%
  filter(inferred_repunit != repunit) %>%
  filter(scaled_likelihood > 0.95) %>%
  filter(!collection %in% c("carnatus", "chrysomelas"))

That made all of the other strong misassignments, with the exception of R013353 go away.

Based on the PCA, let’s take a look at misassigned caurinus samples:

new_self %>%
  filter(collection == "caurinus") %>%
  filter(scaled_likelihood > 0.5) %>%
  arrange(z_score)

There are four caurinus samples with z-scores < -2; one of which is -12! This gives credence to the idea that they are another species, maybe not in the baseline? Or hybrids?

Still not sure what to think about the one that assigned to carnatus. Check meta data?

For the assignment accuracy:

new_self %>%
  filter(repunit == inferred_repunit) %>%
  filter(scaled_likelihood > 0.95)

1579/1647

95.87% accurately assigned at the 95% threshold

new_self %>%
  filter(collection %in% c("carnatus", "chrysomelas")) %>%
  filter(scaled_likelihood > 0.95)

287/346 total GBY

287/346

Lower accuracy for GBY (83%)

Were all mis-id’s to the other species?

new_self %>%
  filter(collection %in% c("carnatus", "chrysomelas")) %>%
  filter(scaled_likelihood > 0.5) %>%
  filter(!inferred_collection %in% c("carnatus", "chrysomelas"))

A single sample assigned to caurinus, but with a z-score of -10!

z-score outliers

Now I’m curious how many of these samples may be from species outside of our baseline according to the z-scores

new_self %>%
  filter(scaled_likelihood > 0.5) %>%
  filter(z_score < -3) %>%
  #filter(z_score > 3) %>%
  arrange(z_score) %>%
  select(indiv, repunit, inferred_repunit, scaled_likelihood, z_score) %>%
  left_join(., meta, by = c("indiv" = "NMFS_DNA_ID")) %>%
  select(-PHENOTYPE, HATCHERY_MARK, TAG_NUMBER, ESTIMATED_DATE, LANDFALL_PORT, CRUISE, HAUL)

Hmmm - a bunch of juveniles. What happens if I remove juvenile samples from the analysis?

For most of these samples, they assign to their reported unit even though they have low z-scores.

Test out making a genepop file without the z-score outliers

Just to test it out, I’m going to remove these outliers from the baseline and then perform rubias’s mixture assignment with them.

Now that the data are separated and formatted, do the mixture assignment

mix_zs

$mixing_proportions

$indiv_posteriors

$mix_prop_traces

$bootstrapped_proportions

NA

mix_zs$indiv_posteriors %>%
  filter(PofZ > 0.5) %>%
  select(mixture_collection, collection, PofZ, z_score, n_miss_loci) %>%
  filter(z_score > -3)

single GBY reporting unit

Since all of the intermediate assignments and 8/9 of the high confidence misassignments are GBY, let’s use a single reporting unit for them.

Now try self-assignment with the single reporting unit

gby_repu_assigned <- gby_repu_2col %>%
  self_assign(., gen_start_col = 5)

Summary Statistics:

1647 Individuals in Sample

90 Loci: Plate_1_A01_Sat_GW603857_consensus.1, Plate_1_A11_Sat_GE820299_consensus.1, Plate_2_A09_Sat_EW986980_consensus.1, Plate_2_C08_Sat_EW987116_consensus.1, Plate_2_G06_Sat_EW987118_consensus.1, Plate_3_C03_Sat_GE798118_consensus.1, Plate_4_E10_Sat_EW976030_consensus.1, Plate_4_G06_Sat_EW976181_consensus.1, tag_id_1049.1, tag_id_108.1, tag_id_1184.1, tag_id_1229.1, tag_id_1272.1, tag_id_1366.1, tag_id_1428.1, tag_id_143.1, tag_id_1441.1, tag_id_1449.1, tag_id_1471.1, tag_id_1498.1, tag_id_1558.1, tag_id_1576.1, tag_id_1598.1, tag_id_1604.1, tag_id_1613.1, tag_id_162.1, tag_id_1652.1, tag_id_170.1, tag_id_1708.1, tag_id_1748.1, tag_id_1751.1, tag_id_1762.1, tag_id_179.1, tag_id_1804.1, tag_id_1808.1, tag_id_1810.1, tag_id_1836.1, tag_id_1850.1, tag_id_1880.1, tag_id_1889.1, tag_id_1915.1, tag_id_1950.1, tag_id_1961.1, tag_id_1966.1, tag_id_1982.1, tag_id_1994.1, tag_id_1999.1, tag_id_2008.1, tag_id_2009.1, tag_id_2017.1, tag_id_2062.1, tag_id_2082.1, tag_id_2114.1, tag_id_2134.1, tag_id_2155.1, tag_id_2178.1, tag_id_2182.1, tag_id_220.1, tag_id_2203.1, tag_id_221.1, tag_id_2214.1, tag_id_2237.1, tag_id_2247.1, tag_id_2258.1, tag_id_2301.1, tag_id_2319.1, tag_id_2368.1, tag_id_2499.1, tag_id_250.1, tag_id_2607.1, tag_id_2635.1, tag_id_265.1, tag_id_325.1, tag_id_402.1, tag_id_410.1, tag_id_436.1, tag_id_55.1, tag_id_572.1, tag_id_67.1, tag_id_770.1, tag_id_788.1, tag_id_843.1, tag_id_855.1, tag_id_874.1, tag_id_875.1, tag_id_879.1, tag_id_913.1, tag_id_942.1, tag_id_981.1, tag_id_987.1

53 Reporting Units: aleutianus, alutus, atrovirens, auriculatus, aurora, babcocki, borealis, chrysomelas, caurinus, chlorostictus, constellatus, crameri, dallii, diaconus, diploproa, elongatus, emphaeus, ensifer, entomelas, flavidus, goodei, hopkinsi, jordani, levis, maliger, melanops, melanostictus, melanostomus, miniatus, moseri, mystinus, nebulosus, nigrocinctus, oculatus, ovalis, paucispinis, pinniger, polyspinis, proriger, rastrelliger, reedi, rosaceus, ruberrimus, rubrivinctus, rufinanus, rufus, saxicola, semicinctus, serranoides, serriceps, umbrosus, wilsoni, zacentrus

54 Collections: aleutianus, alutus, atrovirens, auriculatus, aurora, babcocki, borealis, carnatus, caurinus, chlorostictus, chrysomelas, constellatus, crameri, dallii, diaconus, diploproa, elongatus, emphaeus, ensifer, entomelas, flavidus, goodei, hopkinsi, jordani, levis, maliger, melanops, melanostictus, melanostomus, miniatus, moseri, mystinus, nebulosus, nigrocinctus, oculatus, ovalis, paucispinis, pinniger, polyspinis, proriger, rastrelliger, reedi, rosaceus, ruberrimus, rubrivinctus, rufinanus, rufus, saxicola, semicinctus, serranoides, serriceps, umbrosus, wilsoni, zacentrus

6.6140457397288% of allelic data identified as missing

gby_repu_assigned %>%
  filter(repunit == inferred_repunit) %>%
  #filter(scaled_likelihood > 0.95) %>%
  filter(repunit %in% c("carnatus", "chrysomelas"))

1586/1647

[1] 0.962963

That only bumps it up to 96.3% accurate assignment because there are still lesser assignments (> 0.95)

playing with int hz

int_hz <- t1 %>% 
  group_by(collection, indiv, Locus) %>% 
  summarise(hz = allele[1]!=allele[2]) %>% 
  filter(!is.na(hz)) %>% 
  group_by(indiv, collection) %>% 
  summarise(num_loc = n(), num_hz = sum(hz), fract_hz = num_hz / num_loc)

And now let us plot that.

This is really only useful if the colors of phylogenetically relevant.

Or I could take a subset of the species:

What about making a series of boxplots?

box <- ggplot(int_hz) +
  geom_boxplot(aes(x = collection, y = fract_hz, fill = collection), position = "dodge") +
  theme_bw() +
  xlab("Species") +
  ylab("Internal heterozygosity")
box +
  theme(
    axis.text.x = element_text(angle = 90)
  )

Summary Statistics The lower and upper hinges correspond to the first and third quartiles (the 25th and 75th percentiles).

The upper whisker extends from the hinge to the largest value no further than 1.5 x IQR from the hinge (where IQR is the inter-quartile range, or distance between the first and third quartiles). The lower whisker extends from the hinge to the smallest value at most 1.5 x IQR of the hinge. Data beyond the end of the whiskers are called “outlying” points and are plotted individually.

I still want to set the species phylogenetically rather than alphabetically!

Can we also examine the difference from the mean hz of atrovirens?

int_hz %>%
  group_by(collection) %>%
  summarise(mean_hz = mean(fract_hz)) %>%
  arrange(desc(mean_hz)) %>%
  mutate(diff = mean_hz/0.46478647) %>%
  mutate(perc = diff*100) %>%
  filter(diff < 0.5)